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ABSTRACT 

We compare observations of the high redshift galaxy population to the predictions of the 
galaxy formation model of Croton et al. (2006) and De Lucia & Blaizot (2006). This model, 
implemented on the Millennium Simulation of the concordance ACDM cosmogony, intro- 
duces "radio mode" feedback from the central galaxies of groups and clusters in order to 
obtain quantitative agreement with the luminosity, colour, morphology and clustering prop- 
erties of the present-day galaxy population. Here we construct deep light cone surveys in 
order to compare model predictions to the observed counts and redshift distributions of dis- 
tant galaxies, as well as to their inferred luminosity and mass functions out to redshift 5. With 
the exception of the mass functions, all these properties are sensitive to modelling of dust 
obscuration. A simple but plausible treatment agrees moderately well with most of the data. 
The predicted abundance of relatively massive (^ M*) galaxies appears systematically high at 
high redshift, suggesting that such galaxies assemble earlier in this model than in the real Uni- 
verse. An independent galaxy formation model implemented on the same simulation matches 
the observed mass functions slightly better, so the discrepancy probably reflects incomplete 
galaxy formation physics rather than problems with the underlying cosmogony. 

Key words: galaxies: general - galaxies: formation - galaxies: evolution - galaxies: lumi- 
nosity function, mass function 



1 INTRODUCTION 

J Recent work has used the very large Millennium Simulation to fol- 
■ low the evolution of the galaxy population throughout a large vol- 
ume of the concordance ACDM cosmogony (Springel et al. 2005; 
, Croton et al. 2006; Bower et al. 2006; De Lucia et al. 2006; De Lu- 
' cia & Blaizot 2006). By implementing "semi-analytic" treatments 
of baryonic processes on the stored merger trees of all halos and 
subhalos, the formation and evolution of about 10^ galaxies can be 
simulated in some detail. The inclusion of "radio mode" feedback 
from the central galaxies in groups and clusters allowed these au- 
thors to obtain good fits to the local galaxy population and to cure 
several problems which had plagued earlier galaxy formation mod- 
elling of this type. In particular, they were able to produce galaxy 
luminosity functions with the observed exponential cutoff, domi- 
nated at bright magnitudes by passively evolving, predominantly 
elliptical galaxies. At the same time, this new ingredient provided 
an energetically plausible explanation for the failure of "cooling 
flows" to produce extremely massive galaxies in cluster cores. Most 
of this work compared model predictions to the systematic proper- 
ties and clustering of the observed low redshift galaxy population, 
or studied the predicted formation paths of massive galaxies. Only 
Bower et al. (2006) compared their model in detail to some of the 
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currently available data at high redshift. In the present paper we 
compare these same data and others to the galaxy formation model 
of Croton et al. (2006) as updated by De Lucia & Blaizot (2006) 
and made publicly available through the Millennium Simulation 
data site.^ 

Many recent observational studies have emphasised their de- 
tection of substantial populations of massive galaxies out to at 
least redshift 2 and have seen this as conflicting with expectations 
from hierarchical formation models in the ACDM cosmogony (e.g. 
Cimatti et al. 2002b; Im et al. 2002; Pozzetti et al. 2003; Kashikawa 
et al. 2003; Chen et al. 2003; Somerville et al. 2004). This notion 
reflects in part the fact that early hierarchical models assumed an 
Einstein-de-Sitter cosmogony in which recent evolution is stronger 
than for ACDM (e.g. Fontana et al. 1999), in part an underassess- 
ment of the predictions of the contrasting toy model in which mas- 
sive galaxies assemble at high redshift and thereafter evolve in lu- 
minosity alone (see Kitzbichler & White 2006). Bower et al. (2006) 
find their model to be in good agreement with current observational 
estimates of the abundance of massive galaxies at high redshift, 
while our comparisons below suggest that the model of Croton et al. 
(2006) and De Lucia & Blaizot (2006) appears, if anything, to over- 
predict this abundance. As shown by these authors and particularly 
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by De Lucia et al. (2006) both models predict "anti-hierarchical" 
behaviour, in that star formation completes earlier in more massive 
galaxies. This behaviour clearly does not conflict with the underly- 
ing hierarchical growth of structure in a ACDM cosmogony. 

The current paper is organised as follows. In Section 2 we 
briefly describe the Millennium Simulation and the fiducial galaxy 
formation model we are adopting. Where we have made modifica- 
tions, most significantly in the dust treatment, these are described 
in detail. We also give a detailed account of how we construct mock 
catalogues of galaxies along the backward lightcone of a particu- 
lar simulated field of observation. Many of our methods resemble 
those which Blaizot et al. (2005) implemented in their MOMAF fa- 
cility in order to enable mock observations of simulated galaxy cat- 
alogues of the same type as (though smaller than) the Millennium 
Run catalogues we use here. Our results are summarised in Sec- 
tion 3 where we compare number counts as a function of apparent 
magnitude and redshift with the currently available observational 
data. We also compare the predicted evolution of the luminosity 
and stellar mass functions to results derived from recent observa- 
tional surveys, and we illustrate how the population of galaxies is 
predicted to shift in the colour-absolute magnitude plane. Finally in 
Section 4 we interpret our findings and present our conclusions. 



2 THE MODEL 

2.1 The Millennium dark matter simulation 

We make use of the Millennium Rvm, a very large simulation 
which follows the hierarchical growth of dark matter structures 
from redshift z = 127 to the present. The simulation assumes 
the concordance ACDM cosmology and follows the trajectories of 
2160^ ~ 1.0078 X 10" particles in a periodic box 500Mpc//i on 
a side. A full description is given by Springel et al. (2005); here we 
summarise the main simulation characteristics as follows: 

The adopted cosmological parameter values are consistent 
with a combined analysis of the 2dFGRS (CoUess et al. 2001) 
and the first-year WMAP data (Spergel et al. 2003; Seljak et al. 
2005). Specifically, the simulation takes Qm = f^dm + f^b = 0.25, 
fib = 0.045, h = 0.73, = 0.75, n = 1, and as = 0.9 
where all parameters are defined in the standard way. The adopted 
particle number and simulation volimie imply a particle mass of 
8.6 X 10* /i~^M0. This mass resolution is sufficient to resolve the 
haloes hosting galaxies as faint as 0.1 L* with at least ~ 100 par- 
ticles. The initial conditions nt z — 127 were created by displac- 
ing particles from a homogeneous, 'glass-like' distribution using a 
Gaussian random field with the ACDM linear power spectrum. 

In order to perform such a large simulation on the available 
hardware, a special version of the GADGET-2 code (Springel et al. 
2001b; Springel 2005) was created with very low memory con- 
sumption. The computational algorithm combines a hierarchical 
multipole expansion, or 'tree' method (Barnes & Hut 1986), with 
a Fourier transform particle-mesh method (Hockney & Eastwood 
1981). The short-range gravitational force law is softened on co- 
moving scale 5 /i~^kpc which may be taken as the spatial resolu- 
tion limit of the calculation, thus achieving a dynamic range of 10^ 
in 3D. Data from the simulation were stored at 63 epochs spaced 
approximately logarithmically in time at early times and approxi- 
mately linearly in time at late times (with At ~ SOOMyr). Post- 
processing software identified all resolved dark haloes and their 
subhaloes in each of these outputs and then Unked them together 
between neighboring outputs to construct a detailed formation tree 



for every object present at the final time. Galaxy formation mod- 
elling is then carried out in post-processing on this stored data 
structure. 



2.2 The basic semi-analytic model 

Our semi-analytic model is that of Croton et al. (2006) as updated 
by De Lucia & Blaizot (2006) and made public on the Millenium 
Simulation data download site (see Lemson et al. 2006). These 
models include the physical processes and modelling techniques 
originally introduced by White & Frenk (1991); Kauffmann et al. 
(1993); Kauffmann & Chariot (1998); Kauffmann et al. (1999); 
Kauffmann & Haehnelt (2000); Springel et al. (2001a) and De Lu- 
cia et al. (2004), principally gas cooling, star formation, chemical 
and hydrodynamic feedback from supernovae, stellar population 
synthesis modelling of photometric evolution and growth of super- 
massive black holes by accretion and merging. They also include 
a treatment (based on that of Kravtsov et al. 2004) of the suppres- 
sion of infall onto dwarf galaxies as consequence of reionisation 
heating. More importantly, they include an entirely new treatment 
of "radio mode" feedback from galaxies at the centres of groups 
and clusters containing a static hot gas atmosphere. The equations 
specifying the various aspects of the model and the specific param- 
eter choices made are listed in Croton et al. (2006) and De Lucia & 
Blaizot (2006). The only change made here is in the dust model as 
described in the next section. 



2.3 Improved dust treatment for the fiducial model 

Even at low redshifts, a crucial ingredient in estimating appropriate 

magnitudes for model galaxies, particularly in the B-band, is the 
dust model. For the present-day luminosity function (LF) a simple 
phenomenological treatment calibrated using observations in the 
local universe has traditionally given satisfactory results (Kauff- 
mann et al. 1999). However, the situation at high redshift is more 
delicate because of the much higher predicted gas (and thus dust) 
columns, the highly variable predicted metallicities, and the shorter 
emitted wavelengths corresponding to typical observed photomet- 
ric bands. We found we had to adopt a new approach in order to 
be consistent with current data on extinction in high redshift galax- 
ies. Devriendt et al. (1999) advocate a dust model based on the 
HI column density in the galaxy disk, a quantity that can be esti- 
mated from the cold gas mass and the disk size of a galaxy, both of 
which are available for each galaxy in our semi-analytic model. A 
plausible scaling of dust-to-gas ratio with metallicity can easily be 
incorporated using the metal content given by a chemical evolution 
model (cf. Devriendt & Guiderdoni 2000). Based on this we get 



z 



A 



{Nh) 



2.1 X 102icm-2 



with the average hydrogen column density obtained from 

M 

{Nh) = 



1.4/imp7rr( 



(1) 



(2) 



A/Av here is the extinction curve from Cardelli et al. (1989). We 
assume the dust-to-gas ratio to scale with metallicity and redshift 
as Tjz = {l + z)~i (Zgas/.^0)^ where s = 1.35 for A < 2000 A 
and s = 1.6 for A > 2000 A. The factor oi {1 + zyi in this 
formula is adopted in order to reproduce results for Lyman-break 
galaxies at z ~ 3. Adelberger & Steidel (2000) find (t)i6oo < 2 
at rest-frame 1600 A, showing that dust-to-gas ratios are lower at 
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this redshift compared to the local universe for objects of the same 
Lboi and metallicity (a result echoed in Reddy et al. 2006). This 
behaviour also agrees with a recent study of the dust-to-gas/dust- 
to-metallicity ratio by Inoue (2003). Please note that the average 
extinction of our model galaxies still increases strongly with red- 
shift due both to the ever shorter rest-frame bands we probe and to 
the smaller disk sizes we predict at higher redshift (see equations 1 
and 2 above). 

2.4 Making mock observations: lightcones 

From a theoretical point of view it would be most convenient to 
compare predictions for the basic physical properties of galaxies 
directly with observation, but in practice this is rarely possible. For 
faint and distant object the most observationally accessible proper- 
ties are usually fluxes in specific observer-defined bands. Quantities 
such as stellar mass or star-formation rate (often even redshift) must 
be derived from these quantities and are subject to substantial un- 
certainties stemming primarily from the assumptions on which the 
conversion is based. Moreover which galaxies can be observed at 
all (and so are included in observational samples) is typically con- 
trolled by observational selection effects on apparent magnitude, 
colour, surface brightness, proximity to other images and so on. 

In order to minimise these uncertainties when drawing astro- 
physical conclusions about the galaxy population, it is beneficial 
to have a simulated set of galaxies with known intrinsic properties 
from which "observational" properties can be calculated, and to ap- 
ply the same conversions and selection effects to this mock sample 
as to the real data. One can then assess the accuracy with which the 
underlying physical properties can be inferred. In this approach the 
uncertain relations between fundamental and observable quantities 
become part of the model, and their influence on any conclusions 
drawn can be assessed by varying the corresponding assumptions 
throughout their physically plausible range. A disadvantage is that 
shortcomings in, for example, the galaxy formation model are con- 
volved with many other effects (for example the conversion from 
mass to luminosity) and separation of these effects can be difficult. 
In particular, it may become difficult to identify why a particular 
model disagrees with the data, since effects from many different 
sources may be degenerate. 

We make mock observations of our artificial universe, con- 
structed from the Milleimium Simulation, by positioning a virtual 
observer at zero redshift and finding those galaxies which lie on his 
backward light cone. The backward light cone is defined as the set 
of all light-like worldlines intersecting the position of the observer 
at redshift zero. It is thus a three-dimensional hypersurface in four- 
dimensional space-time satisfying the condition that light emitted 
from every point is received by the observer now. Its space-like pro- 
jection is the volume within the observer's current particle horizon. 
From this sphere, which would correspond to an all-sky observa- 
tion, we cut out a wedge defined by the assumed field-of-view of 
our mock observation. It is common practice to use the term light 
cone for this wedge rather than for the full (all-sky) light cone, and 
we will follow this terminology here. 

The issues which arise in constructing such light cones have 
been addressed in considerable detail by Blaizot et al. (2005). In 
the following we adopt their proposed solutions in some cases (for 
example, when interpolating the photometric properties of galax- 
ies to redshifts for which the data were not stored) and alternative 
solutions in others (for example when dealing with the limitations 
arising from the finite extent of the simulation). We refer readers to 
their paper for further discussion and for illustration of the size of 



the artifacts which can result from the limitations of this construc- 
tion process. 

There are two major problems to address when constructing 
a light cone from the numerical data. The first arises because the 
Milleimium Simulation was carried out in a cubic region of side 
500Mpc//i whereas the comoving distance along the past light 
cone to redshift 1 is 2390Mpc//i and to redshift 6 is 6130Mpc//i. 
Thus deep light cones must use the underlying periodicity and tra- 
verse the fundamental simulation volume a number of times. Care 
is needed to minimise multiple appearances of individual objects, 
and to ensure that when they do occur they are at widely different 
redshifts and are at different positions on the virtual sky. The sec- 
ond problem arises because redshift varies continuously along the 
past light cone whereas we have stored the positions, velocities and 
properties of our galaxies (and of the associated dark matter) only 
at a finite set of redshifts spaced at approximately 300 Myr intervals 
out to z = 1 and progressively closer at higher redshift. We now 
present our adopted solutions to each of these problems in turn. 

2.4.1 How to avoid making a Mleidoscope 

The underlying scale of the Millennium Simulation 500Mpc//i, 
corresponds to the comoving distance to z ~ 0.17. However, we 
want to produce galaxy catalogues which are at least as deep as the 
current observations, and, in practice, to be one or two generations 
in advance. Although the periodicity of the simulation allows us 
to fill space with any required number of replications of the fun- 
damental volume, this leads to obvious artifacts if the simulation 
is viewed along one of its preferred axes. We can avoid this kalei- 
doscopic effect by orienting the survey field appropriately on the 
virtual sky with respect to the three directions defined by the sides 
of the fundamental cube. The "best" choice depends both on the 
shape and depth of the survey being simulated, and on the criteria 
adopted to judge the seriousness of the artifacts to be minimised. 
Here we do not give an optimal solution to the general problem, 
but rather a solution which works acceptably well for deep surveys 
of relatively small fields. 

Consider a cartesian coordinate system with origin at one cor- 
ner of the fundamental cube and with axes parallel to its sides. 
Consider the line-of-sight from this origin passing through the 
point (i/m, L/n, L) where m and n are integers with no com- 
mon factor and L is the side of the cube. This line-of-sight will 
first pass through a periodic image of the origin at the point 
(nL, mL, nmL), i.e. after passing through nm replications of 
the simulation. If we take the observational field to be defined 
by the lines-of sight to the four points ((n ± 0.5/m)L, (m ± 
0.5/n)L, nmL), it will be almost rectangular and it will have total 
volume L'^/S out to distance [v? -\-m? -\--n?m?Y''' L. Furthermore 
no point of the fundamental cube is imaged more than once. This 
geometry thus gives a mock light cone for a near-rectangular survey 
of size x (in radians) with the first duplicate point 

at distance ~ mnL. For example, if we take m = 2 and n = 3 we 
can make a mock light cone for a 4.8° x 3.2° field out to z = 1.37 
without any duplications. For m = 3 and n = 4 we can do the 
same for a 1.6° x 1.2° area out to 2 = 5.6. Choosing m = 1 and 
n = 5 results in a 11.5° x 2.3° survey with no duplications out to 
z = 1.06. 

If we wish to construct a mock survey for a larger field or to a 
greater distance than these numbers allow, then we have to live with 
some replication of structure. Choosing the central line-of-sight of 
to be in a "slanted" direction of the kind just described with m and 
n values matched roughly to the shape of the desired field usually 
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results in large separations of duplicates in angle and/or in redshift. 
Careful optimisation is needed for any specific survey geometry in 
order to get the best possible results. Note that any point within the 
fundamental cube can be chosen as the origin of a mock survey, 
and that, in addition, there are four equivalent central lines-of sight 
around each of the three principal directions of the simulation. It is 
thus possible to make quite a number of equivalent mock surveys 
of a given geometry and so to ensure that the fuU statistical power 
of the Millennium Run is harnessed when estimating statistics from 
these mock surveys. 

Taking into account the above considerations, we select the 
central Une-of sight to be in the direction of the unit vector us de- 
fined by 

/ 2 I 2 , 2 2\l/2 / N 

(m + n +m n ) ' us = (n, m, mn), (3) 

we define a second unit vector ui to be perpendicular both to us 

and to the unit vector along the coordinate direction associated with 
the smaller of m and n (the a;-axis in the above examples) and we 
take a third unit vector U2 to be perpendicular to the first two so 
as to define a right-handed cartesian system. If we define a and 5 
as local angular coordinates on the sky in the directions of ui and 
U2 respectively, with origin in our chosen central direction, then a 
particular 3-dimensional position x corresponds to 

tan a = X • ui /x • ua 
tan 5 = X • U2 /x • ua 

The position x lies within our target rectangular field provided 

I tan a| ^ tan Aa/2 
I tan (5 1 < tanA(5/2, 

Where Aa and AS give desired angular extent of the field in the 
two orthogonal directions (with Aa ^ AS assumed here). Note 
that this formulation of the condition to be within the Ught cone 
does not require any transcendental functions to be applied to the 
galaxy positions, allowing membership to be evaluated efficiently. 
This can be a significant computional advantage when one is re- 
quired to loop over many replications of the (already large) Millen- 
niiun galaxy catalogues. 

We point out in passing that only for comoving coordinates 
within a flat universe do we have the luxury of cutting out light 
cones from our (replicated) simulation volume simply as we would 
in Euclidian geometry. In general this is a much less trivial endeav- 
our that requires accounting for the curvature of the universe as 
well as its expansion with time. (In addition, second-order effects 
like gravitational lensing should, in principle, be taken into account 
for any geometry.) 

2.4.2 How to get seamless transitions between snapshots 

After determining the observer position and survey geometry we fill 
three-dimensional Euclidian space-time with a periodically repli- 
cated grid of simulation boxes, keeping only those which intersect 
our survey. In practice, since the Millennium Simulation data at 
each time are stored in a set of 512 spatially disjoint cells, we keep 
only those cells which intersect the survey. In principle, a galaxy 
within our survey at comoving distance D from the observer should 
be seen as it was at redshift z where 



A problem arises, however, because the positions, velocities and 
physical properties of our galaxies are stored only at a discrete set 



of redshifts Zi corresponding to a discrete set of distances Di. (For 
definiteness we adopt = -Di = and Zi > Zi-i,Di > Dj-i.) 
The comoving distance between outputs is 80 to 240 Mpc/fe, corre- 
sponding to 100 to 380 Myr, depending on redshift. 

One way to deal with this problem would be to interpolate 
the positions, velocities and physical properties of the galaxies at 
each distance D from the output redshifts which bracket it, e.g. 
Zi and Zi+i where Di+i > D > Di. We decided against this 
procedure for several reasons. In the first place, the Millennium 
Simulation appears to give dynamically consistent results for the 
galaxy distribution down to scales of lOkpc or so (see, for exam- 
ple, the 2-point correlation functions in Springel et al. 2006). On 
such scales characteristic orbital timescales are smaller than the 
spacing between our outputs, so interpolation would produce dy- 
namically incorrect velocities and would diffuse structures. In ad- 
dition, the physical properties of the galaxies are not easily inter- 
polated because of impulsive processes such as mergers and star- 
bursts. Rather than interpolating, we have chosen to assign the po- 
sitions, velocities and physical properties stored at redshift Zi to 
all survey galaxies with distances from the observer in the range 
{Di + A+i)/2 > D > {Di + A-i)/2. Individual small scale 
structures are then dynamically consistent throughout this range, 
and the physical properties of the galaxies are offset in time from 
the correct values by at most half of the time spacing between out- 
puts. 

After coarsely filling the volume around the observed light 
cone with simulation cells in this way one can simply chisel off 
the protruding material, i.e. drop all galaxies which do not lie in 
the field according to the condition in Eqn. 4 or which don't satisfy 
{Di + A+i)/2 > D > (Di + A-i)/2. The latter condition 
causes an additional difficulty since galaxies move between snap- 
shots and thus it can happen that a galaxy traverses the imaginary 
boundary {Di + /)f+i ) /2 between the times corresponding to Zi+i 
and Zi. This results in this galaxy being observed either twice or not 
at all, depending on the direction of its motion. We overcome this 
problem for galaxies close to the boundary by linearly interpolat- 
ing their positions between Zi+i and Zi in order to get estimated 
positions at the redshift corresponding to {Di + Di+i)/2. Those 
galaxies whose estimated positions are on the low redshift side of 
the boundary are assigned properties corresponding to Zi, those on 
the high redshift side properties corresponding to Zi+i. 

2.4.3 Getting the right magnitudes 

The observed properties of a galaxy depend not only on its intrin- 
sic physical properties but also on the redshift at which it is ob- 
served. In particular, the apparent magnitudes of galaxies are usu- 
ally measured through a filter with fixed transmission curve in the 
observer's frame. This transmission curve must be blue-shifted to 
each galaxy's redshift and then convolved with the galaxy's spectral 
energy distribution in order to obtain an absolute luminosity which 
can be divided by the square of the luminosity distance to obtain 
the observed flux. A difficulty arises because quantities like abso- 
lute luminosities are accumulated, based on the prior star formation 
history of each object, at the time the semi-analytic simulation is 
carried out, and they are stored in files which give the properties of 
every galaxy at each output redshift Zi. At this stage the light cone 
surveys are not yet defined, so we do not know the redshift at which 
any particular galaxy will be observed in a particular mock survey. 
We are thus unable to define the filter function through which its 
luminosity should be accumulated in order to reproduce properly 
the desired observer-frame band. 
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We deal with this problem in the way suggested by Blaizot 
et al. (2005). We define ahead of time the observer frame magni- 
tudes we wish to predict, for example, Johnson B. When carrying 
out the semi-analytic simulation we then accumulate for all galax- 
ies at redshift Zi not only the absolute magnitude through the B- 
filter blue-shifted to Zi but also those for the same filter shifted to 
the frequency bands corresponding to and Zi+i. For galaxies 
in our mock survey whose physical properties correspond to Zi but 
which appear on the light cone at z > Zi we linearly interpolate an 
estimate for the observer-frame B absolute magnitude (at redshift 
z) between the values stored for filters blue-shifted to Zi and Zi+i. 
Similarly, for those similar galaxies which appear at z; < Zi we 
interpolate the absolute magnitude between the values stored for 
filters blue-shifted to Zi and Zi^i.ll turns out that this interpolation 
is quite important. Without it, discontinuities in density are readily 
apparent in the distribution of simulated galaxies in the observed 
colour-apparent magnitude plane. 

We conclude this chapter by presenting in Fig. 1 an illus- 
trative example, the simulated light cone of a deep survey (to 
Ks(AB) < 24) of a 1.4° x 1.4° field out to z = 3.2. Here intensity 
corresponds to the logarithmic density and the colour encodes the 
offset from the evolving red sequence at the redshift of observation 
(assuming passive evolution after a single burst at z — 6). Large- 
scale structure is evident and is well sampled out to redshifts of at 
least 2 ~ 3 and it is interesting that at z > 2 the reddest galaxies 
are predicted to be in the densest regions even though, as we see be- 
low, many of them are predicted to be dusty strongly star forming 
objects. Individual bright galaxies are predicted to be visible out to 
z ~ 5 in the full light cone. 



3 RESULTS 

In this section we first compare our model to directly measured 
properties of real samples such as their distribution in apparent 
magnitude and redshift. We then consider derived properties which 
require an increasing number of additional assumptions, moving 
from the evolution of rest-frame luminosity functions to that of 
stellar mass distributions. Finally we illustrate the large changes 
predicted for the distribution of galaxies in rest-frame colour and 
absolute magnitude over the redshift range < z < 3. This gives a 
good impression of the interplay between the various mechanisms 
that determine the luminosity and colour of galaxies in our model. 

All magnitudes are in the AB system (rather than Vega) unless 
stated otherwise. 



3.1 Number Counts 

In Fig. 2 we compare predicted galaxy counts obtained from a mock 
survey of a 2 0° area to observational counts from a number of dif- 
ferent surveys. In the BRI bands we use counts over a 0.2 0° area 
in the HDF-N direction by Capak et al. (2004). In the Ks band 
we use both the "wide" area (320 arcmin^ distributed over various 
fields) counts of Kong et al. (2006) and the deeper, but smaller area 
counts in the CDF and HDF-S directions (6 and 7.5 arcmin^ re- 
spectively) by Saracco et al. (2001). It is worth noting that for the 
BRI bands we were able to use the filter transmission curves ap- 
propriate for the Subaru survey, whereas for the K band different 
effective transmission curves apply for the different surveys and 
we have not taken this into account. In order to quantify the ef- 
fect of "cosmic variance" (the fact that large statistical fluctuations 
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Figure 2. Predicted galaxy counts per unit area in four bands compared to a 
survey in the HDF-N direction (0.2 0°) in BRI and to a number of "wide 
field" surveys at Ks (total of 320 arcmin^), as well as to deep observations 
at Ks in the CDF and HDF-S directions (6 and 7.5 arcmin^ respectively). 
The grey shaded error bands show Icr field-to-field variations assuming an 
ai'ea of 100 arcmin^, whereas the hatched en'or bands show the expected 
variations for a smaller square field of area ~ 11 arcmin^. 

are expected in surveys of this size not only from counting statis- 
tics but also from large-scale structure along the line-of-sight) we 
split up our 2 mock survey into 72 fields of size 100 arcmin'^. 
The la scatter among counts in these different areas is shown as 
a grey shaded area surrounding the predicted means for BRI and 
for the brighter K magnitudes. For the fainter K magnitudes we 
split our mock survey into smaller subfields, each with an area of 
~ 11 arcmin^. The la variations among these subfields are shown 
by the hatched band surrounding the predicted K counts at fainter 
magnitudes. Note that this procedure may still somewhat underes- 
timate the cosmic variance since the different subfields are not truly 
independent, but all lie within a single 1.4° x 1.4° mock survey. 

In the light of this limitation, and keeping in mind that our 
dust-model is still rather simple, it is quite surprising to see the 
excellent agreement of the data with our predictions in all three 
optical bands. Agreement at K is less good, and there appears to 
be a significant disrepancy faintward of Kab ~ 21. The model 
predicts almost twice as many galaxies as are observed at Kab '~ 
23, although the agreement is again acceptable at Kab ^ 24.5. 
This disagreement appears well outside the statistical errors, but it 
should be borne in mind that K magnitudes are extremely difficult 
to measure at such faint levels, and it is possible that the measured 
quantity does not correspond to the total magnitude assumed in our 
modelling. 

3.2 Redshift Distributions for If -selected samples 

In Fig. 3 we give the redshift distributions predicted for appar- 
ent magnitude limited galaxy samples complete for AT ^ 21.8, 
K ^ 23.3, and K ^ 25.8. We compare the first of these to data 
for a 52 arcmin^ field from K20 (Cimatti et al. 2002b) and for 
a 160 arcmin^ overlapping field from GOODS (Mobasher et al. 
2004). (Note that the name K20 comes from the survey limit in 
the Vega system. The two systems are approximately related by 
Kab ~ -K^Voga + 1.83.). At the intermediate depth we compare 
to the photometric redshift distribution obtained by Caputi et al. 
(2006) for a 131 arcmin^ field in the direction of the Chandra Deep 
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Figure 1. Light cone for a 1.4° X 1.4° field out to z = 3.2. All galaxies above an apparent magnitude limit Ks{AB) < 24 are shown, where intensity 
corresponds to the logarithmic density and the color denotes the offset from the evolving red sequence. 



Field South (CDF-S). For the faintest magnitude limit we compare 
to the photo-z distribution of a much smaller 4 arcmin'^ area in the 
Subaru Deep Field, as obtained by Kashikawa et al. (2003). Again 
we split up our simulated field into sub-fields of size 100 arcmin'^ 
(4 arcmin^ for K ^5 24) in order to get an estimate of the expected 
Icr scatter, which we indicate by grey shaded areas. For these small 
fields cosmic variance is quite substantial and the counts can be in- 
fluenced significantly by individual galaxy clusters. This effect is 
clearly visible in the K20 and GOODS data, where a pronounced 
spike is present at z ~ 0.7. In addition, systematic problems with 
the photometric redshift determinations might distort the redshift 
distributions in some ranges. 

Despite these uncertainties, our model predictions appear 
somewhat high over the redshift range 0.5 < z < 1.5 for the 
K20 and GOODS samples. The deeper K ^ 23.3 observations are 
overpredicted by a factor of 2 to 3 over the range 1 < z < 3. For 
the faintest sample there is an apparent overprediction by a some- 
what smaller factor over this same redshift range. Comparing with 
Fig. 2, we see that the total overprediction at each of these magni- 
tudes is consistent with that seen in the counts themselves, although 
it should be borne in mind that the CDF-S field is common to both 
datasets. The differences we find are larger than the predicted cos- 
mic variance so they presumably indicate problems with the model 



(incorrect physics?) with the observational data (systematics in the 
magnitudes or photo-z's) or both. 

3.3 Luminosity Function evolution 

Croton et al. (2006) demonstrated that at z = the luminosity 
function (LF) for our model agrees well with observation both in 
b,j and in K. Splitting galaxies according to their intrinsic colours, 
these authors also found quite good fits to the LF's for red and blue 
galaxies separately, with some discrepancies for faint red galaxies. 
Here we compare the evolution of the LF predicted by our model 
in rest-frame B and K band with recent observational results. 

3.3. 1 The B-hand Luminosity Function 

In Fig. 4 we compare the evolution of the rest frame _B-band LF 
predicted by our simulation to results from the DEEP2 survey 
(Willmer et al. 2005). As a z = standard we use the local LF 
from the 2dF survey Norberg et al. (2002). This is compared with 
our model in the top-left panel and is repeated as a thin red line 
in each of the other panels, where the high redshift data are indi- 
cated by points with error bars. Our predicted LF is shown in each 
panel as a solid line with a grey area indicating the la scatter to be 
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Figure 3. Predicted redshift distributions for galaxies to a magnitude limit 
of K 21.8 (solid), 23.3 (dotted), and 25.8 (dashed). These are compared 
to observational results from K20 and GOODS (K ^ 21.8), from CDF-S 
{K ^ 23.3) and from SDF (K 25.8). The latter two are derived purely 
from photometric redshift estimates. Error bars on the observational points 
are based on counting statistics only. Grey shaded areas indicate the la 
field-to-field scatter assuming an area of 100 arcmin^ for the two brighter 
magnitude limits and 4 arcmin^ for K ^ 25.8. 

expected for an estimate from a survey similar in effective volume 
to the observational survey. (Note that in all cases the Millennium 
Simulation is much larger than this effective volume, so that count- 
ing noise uncertainties in the prediction are negligible.) 

At z = the agreement between model and observation is 
excellent. This is a consequence of the fact that Croton et al. (2006) 
and De Lucia et al. (2006) adjusted model parameters in order to 
optimise this agreement. However, over the full redshift range from 
z = 0.2 to 1.2 the predicted LP's agree with the DEEP2 data at the 
la level or better. On closer examination, it appears that the model 
somewhat overpredicts the observational abundance fainter than the 
knee of the luminosity function, by a factor ~ 1.5 depending on 
redshift. On the other hand, at the higher redshifts very luminous 
galaxies appear slightly more abundant in the real data than in the 
model. It is important to keep in mind that our dust model has a 
strong influence here, and plausible modifications to it might ac- 
count for either or both of these minor discrepancies. In general, 
the agreement with the data seems quite impressive, at least in this 
band and over this redshift range. 



3.3.2 The K-hand Luminosity Function 

Model predictions for the rest-frame iC-band LF should, in prin- 
ciple, be more robust than predictions for the rest-frame B-band, 
because the effects of our uncertain dust modelling are then much 
weaker. On the other hand, observational determinations of the LF 
at rest-frame K are more uncertain than at rest-frame B, because 
the magnitudes of high redshift galaxies must then be inferred 
by extrapolation beyond the wavelength region directly measured, 
rather than interpolated between the observed bands. This situation 
is improving rapidly as deep data at wavelengths beyond 2fj, be- 
come available from Spitzer. 

As can be seen in Fig. 5, our predictions for the evolution of 
the rest-frame AT-band LF show the same behaviour as for the B- 
band. The local result from Cole et al. (2001) is reproduced well, 
as illustrated in the upper left panel and already demonstrated in 




Figure 4. Comparison of the rest-frame B-band LF predicted by our sim- 
ulation for redshifts in the range < z < 1.2 to observational estimates 
from Norberg et al. (2002) at z = and from Willmer et al. (2005) at higher 
redshifts. The local LF of the upper left panel is repeated as a thin red line 
in each of the other panels. A grey shaded region surrounding each model 
prediction shows the Icr scatter expected in observational estimates based 
on samples similar in size to the corresponding observational sample. 

Croton et al. (2006). The observed z = function is reproduced 
as a thin red line in the other panels in order to make the amount 
of evolution more apparent. At higher redshifts we compare with 
observational determinations from Pozzetti et al. (2003) for the 52 
arcmin^ of the K20 survey, from Feulner et al. (2003) for the 600 
armin^ of the MUNICS sample and from Saracco et al. (2006) for a 
5.5 arcmin^ area in the HDF-S. In these plots we give error bars as 
quoted by the original papers, but we note that these are based on 
counting statistics only and additional uncertainties are expected 
due to clustering, particularly for the smaller fields. Furthermore, 
photo-z's are used for a significant number of galaxies in these de- 
terminations which may lead to additional systematic uncertainties 
in the results. Given the scatter between the various observational 
determinations, the disagreements between model and data do not 
look particularly serious. The models do appear to overpredict the 
abundance of galaxies near the knee of the luminosity function, 
perhaps by a factor of 2 at the highest redshift, echoing the discrep- 
ancies found above when comparing with i^-band galaxy counts 
and redshift distributions. 



3.3.3 Evolution of Luminosity Function parameters 

In order to display the evolution of the luminosity function in our 
models more effectively, we have fit Schechter (1976) functions 
to the simulation data for the rest-frame B and AT-bands at every 
stored output time. In most cases these functions are a good enough 
fit to give a fair representation of the numerical results. In Fig. 6 
we plot the evolution with redshift of the parameters <1>* and Al* 
and of the volume luminosity density, j = ^* L*r{a + 2) using 
thick solid lines, and we compare with fits to observational data. 
For each observational point in the if-band panels we indicate the 
(often broad) redshift range to which it refers by a horizontal bar. 
The vertical bar indicates the uncertainty quoted by the original 
authors. 

Not surprisingly, the results of the last section are confirmed. 
For the model $* increases slightly with redshift out to about 
z = 1.5, whereas the observations imply a relatively steep de- 
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Figure 5. Comparison of the evolution predicted for the rest-frame K-band 
LF to observational determinations from Cole et al. (2001) at low z (upper 
left panel, repeated as a thin red line in the other panels) and from Pozzetti 
et al. (2003), Feulner et al. (2003) and Saracco et al. (2006) at higher red- 
shifts. 



cline over this same redshift range. This holds for both photomet- 
ric bands. For A/* we see brightening both in the models and in 
the observations, but the effect is more pronounced in the latter. 
In K the models predict Al* to be almost independent of redshift. 
Derivations of <I>* and M* from observational data using maxi- 
mum likelihood techniques usually give results where the errors in 
the two quantities correlate in a direction almost parallel to lines of 
constant luminosity density. For this reason we expect j to be more 
robustly determined from the data than either $* or M* individu- 
ally. It is interesting that the apparent deviations between data and 
model for $* and M* largely compensate, so that the model pre- 
dicts an evolution of j which is quite similar to that inferred from 
the observations. This is particularly striking at rest-frame B. At 
rest-frame K the observational error bars are still too large to draw 
firm conclusions, but a non-evolving luminosity density represents 
the data somewhat better than does our model, again confirming the 
conclusions we drew in earlier sections. 



3.4 The evolution of the stellar mass function 

The evolution of the abundance of galaxies as a function of their 
stellar mass is one of the most direct predictions of galaxy for- 
mation models. It depends on the treatment of gas cooling, star- 
formation and feedback, but not directly on the luminous properties 
of the stars or on the dust modelling. (There remains an indirect 
dependence on the latter since observations of galaxy luminosities 
are typically used to set uncertain efficiency parameters in the mod- 
elling.) The stellar masses of galaxies can also be inferred relatively 
robustly from observational data provided sufficient observational 
information is available (e.g. Bell & de Jong 2001, Kauffmann et 
al 2003). At high redshift, however, such inferences become very 
uncertain unless data at wavelengths beyond 2/i are available (e.g. 
from Spitzer). The observationally inferred masses also depend sys- 
tematically on the assumed Initial Mass Function for star formation 
(usually taken to be universal) and it is important to ensure that 
consistent assumptions about the IMF are made when comparing 
observation and theory. 

Bearing in mind these caveats, Fig. 7 compares the mass func- 
tions predicted by our model to local data from Cole et al. (2001) 
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Figure 6. The evolution of the $* and M* parameters of Schechter (1976) 
fits to luminosity functions in rest-frame B and K (with a held constant 
at the values indicated). We also show the evolution of the total luminosity 
density j inferred from these parameters. In each panel the solid line de- 
notes the model prediction and the symbols are data from different sources 
(with potentially different a). The B band data comprises observations 
from Poll et al. (2003, + symbols) and Faber et al. (2005, diamonds), who 
also provide a compilation from the literature (filled circles), whereas the K 
band data are observations (diamonds) and a literature compilation (filled 
circles) from Saracco et al. (2006). 



as well as to high-redshift estimates from Drory et al. (2005), based 
on the MUNICS survey, and from Fontana et al. (2006) based on 
the MUSIC-GOODS data. The latter study uses data in the 3.6 to 
8/1 bands from Spitzer to constrain the spectral energy distributions 
of the galaxies and so should give substantially more reliable re- 
sults at high redshift than the former. In Fig. 7 the model mass 
functions at z > are shown both before and after convolution 
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Figure 7. Evolution of the stellar mass function in the redshift range 
z = — 4.5. Local data are from Cole et al. (2001) and are repeated as a 
black dashed line in the higher redshift panels. High redshift data are taken 
from Drory et al. (2005, symbols) and Fontana et al. (2006, grey shaded 
areas). Model predictions are shown both with (red) and without (black) 
convolution with a normal distribution of standard deviation 0.25 repre- 
senting measurement errors in log A/, . At z = we consider the mass 
determinations precise enough to neglect this effect. 



with a gaussian in log M* with standard deviation 0.25. This is in- 
tended to represent the uncertainty in the observational determina- 
tions of stellar mass. This error may be appropriate for the MUSIC- 
GOODS sample at all redshifts, but it is certainly too small to rep- 
resent uncertainties in the MUNICS mass estimates at high red- 
shift. We note that such errors weaken the apparent strength of the 
quasi-exponential cut-off at high masses. We neglect their effects 
at 2: = 0. 

Our model is nicely consistent with the observed mass func- 
tion in the local Universe, but it clearly overpredicts the abundance 
of galaxies at redshifts between 1 and 3. The observed evolution 
relative to the z = function (indicated by the dashed line in each 
of the higher redshift panels) is strong, while the model prediction 
is rather more modest. In the stellar mass range 10^° to 10" A/0 
where the observational estimats appear most reliable the overpre- 
diction reaches a factor of about 2 at z = 2. This is nicely consis- 
tent with the conclusions we reached in earlier sections based on 
number counts, redshift distributions and luminosity functions, but 
unfortunately the scatter between the various observational deter- 
minations is large enough to prevent any firm conclusion. 



3.5 The evolution of the colour-magnitude relation 

Recent studies of the high redshift galaxy population have often 
stressed the presence of massive objects with colours similar to 
those expected for fully formed and passively evolving ellipticals 
(e.g. Renzini 2006, and references therein). This is usually pre- 
sented as a potential problem for "hierarchical" models of galaxy 
formation where star formation and merging continue to play a ma- 
jor role in the build up of galaxies even at recent times. In order 
to illustrate how these processes are reflected in the colours and 
magnitudes of galaxies in our simulation, we show in Fig. 8 the 
colour-magnitude diagram for 10000 galaxies randomly sampled 
from a 2.5 x 10^/i~^Mpc^ volume at redshifts 2; = 0, 1, 2 and 
3. At z = the well-known bi-modal distribution of colours is 
very evident. A tight red-sequence of passively evolving objects is 
present with a slope reflecting a relation between mass and metal- 
licity. There is also a "blue cloud" of star-forming systems. A suc- 
cess of the model emphasised by Croton et al. (2006) is the fact that 
the brightest galaxies all lie on the red sequence at z — 0. This is 
a consequence of including a treatment of "radio feedback" from 
AGN. 

To allow better appreciation of the evolution to high red- 
shift, we also show logarithmically spaced contours of the colour- 
magnitude distribution of all galaxies in a 1.5 x lO^/i^'^Mpc"^ vol- 
ume as black contours in the panels of Fig. 8. The bluing of the up- 
per envelope with increasing redshift is very clear and is consistent 
with passive evolution of the red sequence. We illustrate this by fit- 
ting a population synthesis model to the ridge line of the z = red 
sequence, assuming a single burst of star formation at 2; = 6 and 
a metallicity which varies with stellar mass. This model is shown 
as a red line not only at 2: = (where it was fit) but also at the 
earlier redshifts. Notice that although there are galaxies with red 
sequence colours at all redshifts, the sequence becomes less and 
less well-defined at earlier times, with a substantial number of ob- 
jects appearing redder than the passively evolving systems. These 
are compact, gas- and metal-rich galaxies where our model predicts 
very substantial amounts of reddening. Recent surveys of distant 
Extremly Red Objects have found substantial numbers of such sys- 
tems (Cimatti et al. 2002a, 2003; Le Fevre et al. 2005; Kong et al. 
2006), but it remains to be seen if our model can account quantita- 
tively for their properties. 

The colours of galaxies in the blue cloud also become bluer 
at high redshift. This is a consequence of an increase in the typical 
ratio of current to past average star formation rate in these galax- 
ies. The difference between star-forming systems and "true" red- 
sequence galaxies becomes blurred at high redshift in our model 
because of the increasingly important effects of dust. 



4 DISCUSSION AND CONCLUSIONS 

The model we have used in this paper is that of Springel et al. 
(2005) and Croton et al. (2006) and as updated by De Lucia & 
Blaizot (2006) and made public through the Millennium Simula- 
tion download site (see Lemson et al. 2006). Earlier work has com- 
pared this model to a wide range of properties of low redshift galax- 
ies: their luminosity functions, their bi-modal luminosity-colour- 
morphology distribution and their TuUy-Fisher relation (Croton 
et al. 2006); their spatial clustering as inferred from two-point cor- 
relations (Springel et al. 2005; Li et al. 2006; Meyer et al. 2006) and 
from fits to halo occupation distribution models (Wang et al. 2006; 
Weinmann et al. 2006); their HI gas content (Meyer et al. 2006); 
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Figure 8. Evolution of the colour-magnitude distribution of galaxies in rest- 
frame B and K over the redshift range z = [0,3]. Randomly selected 
10000 galaxies in a 2.5 X 10^/i~'^Mpc'' volume of the simulation are plot- 
ted in each panel. For comparison purposes the distribution of a much larger 
sample in a 1.5 X lO^/i^^Mpc"^ volume is indicated as logarithmically 
spaced contours. The red solid line indicates the colour magnitude relation 
predicted for stellar populations formed in a single burst at z = 6 and 
evolving passively thereafter. The metallicity of the populations has been 
adjusted as a function of stellar mass to fit the ridge line of the z = red 
sequence. 

and their assembly histories within clusters (De Lucia et al. 2006; 
De Lucia & Blaizot 2006). Although Croton et al. (2006) compared 
the evolution of the global star formation rate and the global black 
hole accretion rate of the model to observation, the current paper 
is the first to compare its predictions in detail with observations of 
high redshift galaxies. 

Our comparison to galaxy counts, to redshift distributions and 
to observational estimates of luminosity and mass functions at high 
redshift paints a consistent picture despite large statistical uncer- 
tainties and some significant technical issues. Our model appears 
to have too many relatively massive galaxies at high redshift and 
these galaxies appear to be too red. Thus, while we fit optical 
galaxy counts well up to densities of 30 gal/mag/arcmin^, we start 
to overpredict numbers in the K band at densities above about 3 
gal/mag/arcmin^. This overabundance of apparently red galaxies 
shows up in the redshift distributions as an overprediction of the 
number of galaxies with A' ~ 23 to 25 at redshift between about 

1 and 3. These correspond to moderately massive systems near the 
knee of the luminosity function, and indeed, while our rest-frame 
B luminosity functions appear compatible with observation out to 
z ~ 1, at rest-frame K our luminosity functions are noticeably 
high beyond z — 0.5 except possibly for the brightest objects. The 
problem shows up most clearly in our mass functions which over- 
predict observationally estimated abundances by about a factor of 

2 at z = 2. Apparently the mass function of galaxies evolved more 
strongly in the real Universe than in our simulation. 

A galaxy formation model with similar basic ingredients to 
ours, but with important differences of detail has been indepen- 
dently implemented on the Millennium Simulation by Bower et al. 
(2006). This model is also publicly available at the download site. It 
fits low redshift galaxy luminosity functions as well as our model, 
but the comparisons which Bower et al. (2006) show to high- 
redshift luminosity and mass function data (essentially the same 
datasets we use here) demonstrate somewhat better agreement than 



we find in this paper. In the mass range 10^" < h^M,/MQ < 10" 
the abundances predicted by their model are lower than ours by 
about 20% at 2 = 1 and by about 30% at both 2 = 2 and 2 = 3.5, 
despite the fact that at 2 = the two models agree very well. This is 
consistent with the fact that their model forms 20% of all its stars by 
2 = 3.2 and 50% by 2 — 1.65 whereas the corresponding redshifts 
for our model are 2 = 3.6 and 2 — 1.9. These differences arise 
from details of the star formation and feedback models adopted in 
the two cases. 

In summary, both the Bower et al. (2006) simulation and our 
own are consistent with most current faint galaxy data. Thus there 
seems no difficulty in reconciling the observed properties of dis- 
tant objects with hierarchical galaxy formation. The fact that pre- 
dictions from the two simulations differ at a level which can be 
marginally separated by the observations, means that currently ac- 
cessible properties of distant galaxies can significantly constrain 
models of this type, and hence the detailed physics which controls 
the formation and the observable properties of galaxies. The fact 
that the model we test here apparently overpredicts the abundance 
of moderately massive galaxies at high redshift, despite the fact 
that late merging plays a major role in the build-up of its more 
massive galaxies (e.g. De Lucia et al. 2006; De Lucia & Blaizot 
2006), demonstrates that current data are still far from constrain- 
ing the importance of this process. As the data improve, the models 
will have to improve also to remain consistent with them. This in- 
terplay between theory and observation should eventually lead to a 
more convincing and more complete picture of how galaxies came 
to take their present forms. 
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